source('00_util_scripts/mod_bplot.R')
library(readxl)

proj.nm <- 'mission/NMOSD-SNP/'

# odds ratios --------
tjsnp <- read_delim('00_util_scripts/data/nmosd.tianjin.2505.txt') |>
  mutate(cohort = str_replace(cohort, 'ianjin', 'ianjin (n=89)') |>
           str_replace('jing', 'jing (n=103)'))

tjg2r <- tjsnp |>
  filter(site == 'G396R')

g2r_1018 <- read_csv('mission/NMOSD-SNP/tianjin-nmosd-1018.csv')

tjg2r <- g2r_1018 |>
  count(Call) |>
  mutate(count = n - 1, genotype = Call, cohort = 'NMOSD Tianjin (n=89)',
         count = ifelse(Call == 'GG', count - 3, count),
         .keep = 'none') |>
  full_join(tjg2r) |>
  summarise(count = sum(count), .by = c(genotype, cohort)) |>
  mutate(cohort = ifelse(str_detect(cohort, 'NMO'), 'NMOSD Tianjin (n=152)',
                         cohort))

tjg2r |>
  ggplot(aes(count, cohort, fill = genotype)) +
  geom_col(position = 'fill') +
  scale_fill_manual(values = c('grey','blue','red')) +
  labs(x = 'Proportion', title = 'IGHG1-G396R genotype: NMOSD vs HC cohort') +
  theme_pubr()

g1 <- last_plot()

tjsnp |>
  filter(site == 'I232T') |>
  ggplot(aes(count, cohort, fill = genotype)) +
  geom_col(position = 'fill') +
  scale_fill_manual(values = c('grey','blue','red')) +
  labs(x = 'Proportion', title = 'FCGR2B-I232T genotype: NMOSD vs HC cohort') +
  theme_pubr()

g2 <- last_plot()

g1/g2

tjg2r |>
  model_OR_forests(ref_group = '1KGP Beijing (n=103)', alt_allele = 'R') |>
  ggplot(aes(or_value, model)) +
  geom_point(size = 6) +
  geom_text(aes(label = str_c('p=',signif(p.value, 2))), nudge_y = .2) +
  geom_linerange(aes(xmin = ci_lower, xmax = ci_upper)) +
  geom_vline(xintercept = 1, linetype = 'dashed') +
  scale_x_log10() +
  theme_pubr() +
  labs(x = 'Odds ratio (95% CI)', subtitle = '152 NMOSD vs 103 HC',
       title = 'IGHG1-G396R correlation with NMOSD')

tjg2r |>
  model_OR_forests(ref_group = '1KGP Beijing (n=103)', alt_allele = 'R')

tjsnp |>
  filter(site == 'I232T') |>
  model_OR_forests(ref_group = '1KGP Beijing (n=103)', alt_allele = 'T') |>
  ggplot(aes(or_value, model)) +
  geom_point(size = 6) +
  geom_text(aes(label = str_c('p=',signif(p.value, 3))), nudge_y = .2) +
  geom_linerange(aes(xmin = ci_lower, xmax = ci_upper)) +
  geom_vline(xintercept = 1, linetype = 'dashed') +
  scale_x_log10(breaks = c(.4,.7,1,1.3,1.6)) +
  theme_pubr() +
  labs(x = 'Odds ratio (95% CI)',
       title = 'FCGR2B-I232T correlation with NMOSD')

## with sex as co-variate --------
i2t_genotype <- query_snp_maf('rs1050501', type = 'genotypes')
g2r_genotype <- query_snp_maf('rs117518546', type = 'genotypes')

### 1KGP pop meta ------------
pop_1kgp <- read_csv('00_util_scripts/ref/1kgp.2504sample.csv')

pop_1kgp |>
  filter(str_detect(sample, '[A-Z]+\\d+'))

pop_1kgp <- i2t_genotype |>
  filter(str_detect(sample, '1000G')) |>
  mutate(sample = str_extract(sample, '[A-Z]+\\d+'),
         gender = NULL, FCGR2B_I232T = genotype, .keep = 'none') |>
  right_join(pop_1kgp)

pop_1kgp <- g2r_genotype |>
  filter(str_detect(sample, '1000G')) |>
  mutate(sample = str_extract(sample, '[A-Z]+\\d+'),
         gender = NULL, IGHG1_G396R = genotype, .keep = 'none') |>
  right_join(pop_1kgp) |>
  write_csv('00_util_scripts/ref/1kgp.2504sample.csv')

chb_1kgp <- pop_1kgp |>
  filter(subpop == 'CHB') |>
  mutate(cohort = '1KGP Beijing')

tjcvr <- read_excel('tianjin-NMOSD-SNP.xlsx')

chb_1kgp |>
  mutate(g2r_ac = str_count(IGHG1_G396R, 'T')) |>
  ggplot(aes(g2r_ac)) + geom_histogram()

cvr2c <- tjcvr |>
  mutate(IGHG1_G396R = `IGHG1-G396R`, FCGR2B_I232T = `FCGR2B-I232T`,
         sex = ifelse(`性别` == '男', 'M', 'F'), cohort = 'NMOSD Tianjin',
         .keep = 'none') |>
  bind_rows(chb_1kgp)

### G2R ---------
cvr2c_g2r <- cvr2c |>
  mutate(cohort, 
         genotype = strrep('R',str_count(IGHG1_G396R, 'R|T')),
         sex = sex == 'M', .keep = 'none') |>
  summarise(count = n(), .by = 1:3)

cvr2c_g2r_or <- cvr2c_g2r |>
  model_OR_forests(ref_group = '1KGP Beijing', alt_allele = 'R')

cvr2c_g2r_sex_adj <- cvr2c_g2r |>
  model_OR_forests(ref_group = '1KGP Beijing', alt_allele = 'R',
                   covariate = '+sex')

cvr2c_g2r_res <- list('None' = cvr2c_g2r_or,
                      'Sex' = cvr2c_g2r_sex_adj) |>
  bind_rows(.id = 'covariate')

cvr2c_g2r_res |>
  write_csv('NMOSD_G2R_OR_250603.csv')

cvr2c_g2r_res |>
  ggplot(aes(or_value, model)) +
  geom_point(size = 6) +
  geom_text(aes(label = str_c('p=',signif(p.value, 3))), nudge_y = .35) +
  geom_linerange(aes(xmin = ci_lower, xmax = ci_upper)) +
  geom_vline(xintercept = 1, linetype = 'dashed') +
  scale_x_log10() +
  theme_pubr() +
  facet_wrap(~covariate, labeller = 'label_both', ncol = 1) +
  labs(x = 'Odds ratio (95% CI)',
       title = 'IGHG1-G396R correlation with NMOSD')

### I2T ---------
cvr2c_i2t <- cvr2c |>
  mutate(cohort, 
         genotype = case_match(cohort,
                               '1KGP Beijing' ~ strrep('T',str_count(FCGR2B_I232T, 'C')),
                               .default = FCGR2B_I232T),
         sex = sex == 'M', .keep = 'none') |>
  summarise(count = n(), .by = 1:3)

cvr2c_i2t_or <- cvr2c_i2t |>
  model_OR_forests(ref_group = '1KGP Beijing', alt_allele = 'T')

cvr2c_i2t_sex_adj <- cvr2c_i2t |>
  model_OR_forests(ref_group = '1KGP Beijing', alt_allele = 'T',
                   covariate = '+sex')

cvr2c_i2t_res <- list('None' = cvr2c_i2t_or,
                      'Sex' = cvr2c_i2t_sex_adj) |>
  bind_rows(.id = 'covariate')

cvr2c_i2t_res |>
  write_csv('NMOSD_I2T_OR_250603.csv')

cvr2c_i2t_res |>
  ggplot(aes(or_value, model)) +
  geom_point(size = 6) +
  geom_text(aes(label = str_c('p=',signif(p.value, 3))), nudge_y = .35) +
  geom_linerange(aes(xmin = ci_lower, xmax = ci_upper)) +
  geom_vline(xintercept = 1, linetype = 'dashed') +
  scale_x_log10() +
  theme_pubr() +
  facet_wrap(~covariate, labeller = 'label_both', ncol = 1) +
  labs(x = 'Odds ratio (95% CI)',
       title = 'FCGR2B-I232T correlation with NMOSD')

## rs146563847 ---------
rs146 <- read_tsv('mission/NMOSD-SNP/rs146563847.tsv')

rs146 |>
  ggplot(aes(`465-510`, `533-580`, color = as.character(Genotype))) +
  geom_point()

rs146_1kgp <- query_snp_maf('rs146563847', type = 'genotypes')

rs146_1kgp |>
  mutate(sample = str_remove(sample, '.+:')) |>
  left_join(pop_1kgp) |>
  filter(subpop == 'CHB') |>
  mutate(genotype = case_match(genotype,
                               'T|T' ~ 'TT',
                               'C|C' ~ 'CC',
                               .default = 'CT'),
         IGHG1_G396R = case_match(IGHG1_G396R,
                                  'T|T' ~ 'TT',
                                  'C|C' ~ 'CC',
                                  .default = 'CT')) |>
  ggplot(aes(genotype, fill = genotype)) +
  geom_bar() +
  geom_text(aes(label = after_stat(count)), stat = 'count', nudge_y = 1) +
  facet_grid(~IGHG1_G396R, labeller = 'label_both') +
  labs(x = 'rs146563847',
       title = 'Linkage between IGHG1-G396R and rs146563847 in 1KGP population') +
  theme_bw()

rs146_chb <- rs146_1kgp |>
  mutate(sample = str_remove(sample, '.+:')) |>
  left_join(pop_1kgp) |>
  filter(subpop == 'CHB') |>
  mutate(genotype = str_count(genotype, 'C'), sex,
         cohort = '1KGP Beijing (n=103)', .keep = 'none') 

rs146_chb

tjcvr <- read_excel('mission/NMOSD-SNP/tianjin-NMOSD-SNP.xlsx')

tjcvr <- tjcvr |>
  mutate(sex = ifelse(性别 == '女', 'F', 'M'), ID, `IGHG1-G396R`, `FCGR2B-I232T`,
         .keep = 'none') |>
  write_csv('mission/NMOSD-SNP/tianjin-89NMOSD-meta.csv')

rs146_2coh <- tjcvr |>
  inner_join(rs146) |>
  mutate(genotype = as.integer(Genotype), cohort = 'NMOSD Tianjin (n=89)', sex,
         .keep = 'none') |>
  bind_rows(rs146_chb)

rs146_2coh |>
  mutate(genotype = case_match(genotype, 0L ~ 'TT',
                               1L ~ 'CT', .default = 'CC') |>
           fct_relevel('TT','CT')) |>
  ggplot(aes(y = cohort, fill = genotype)) +
  geom_bar(position = 'fill') +
  theme_pubr() +
  scale_fill_manual(values = c('grey','royalblue','tomato')) +
  labs(x = 'Proportion', title = 'rs146563847 genotype: NMOSD vs HC cohort')

ggsave('mission/NMOSD-SNP/figs/rs146564847_barplot.png', width = 8, height = 3)

### OR plot ----------
rs146_2coh <- rs146_2coh |>
  count(genotype, cohort, sex) |>
  mutate(genotype = strrep('C', genotype), sex = sex == 'M',
         count = n) 

rs146_or <- rs146_2coh |>
  model_OR_forests(ref_group = '1KGP Beijing (n=103)', alt_allele = 'C')

rs146_sexcor <- rs146_2coh |>
  model_OR_forests(ref_group = '1KGP Beijing (n=103)', alt_allele = 'C',
                   covariate = '+sex')

rs146_2or <- list('None' = rs146_or,
                      'Sex' = rs146_sexcor) |>
  bind_rows(.id = 'covariate')

rs146_2or |>
  write_source_csv('rs146563847_odds_ratio', project = 'mission/NMOSD-SNP/')

rs146_2or |>
  ggplot(aes(or_value, model)) +
  geom_point(size = 6) +
  geom_text(aes(label = str_c('p=',signif(p.value, 3))), nudge_y = .35) +
  geom_linerange(aes(xmin = ci_lower, xmax = ci_upper)) +
  geom_vline(xintercept = 1, linetype = 'dashed') +
  scale_x_log10() +
  theme_pubr() +
  facet_wrap(~covariate, labeller = 'label_both', ncol = 1) +
  labs(x = 'Odds ratio (95% CI)',
       title = 'rs146563847 correlation with NMOSD')

# rs188721059 -----------
tjcvr <-
  read_csv('mission/NMOSD-SNP/tianjin-89NMOSD-meta.csv')

rs188 <- read_delim('mission/NMOSD-SNP/rs188721059.tsv') |>
  mutate(Genotype = case_match(Genotype, '0？1？' ~ '1',
                               'Unknown' ~ '2',
                               .default = Genotype))

rs188 |>
  ggplot(aes(`465-510`, `533-580`, color = Genotype)) +
  geom_point()

rs188_1kgp <- query_snp_maf('rs188721059', type = 'genotypes') |>
  left_join(pop_1kgp) |>
  filter(subpop == 'CHB') |>
  mutate(sex, genotype, cohort = '1KGP Beijing (n=103)', .keep = 'none')

rs188_1kgp

rs188_2coh <- rs188 |>
  inner_join(tjcvr) |>
  mutate(genotype = as.integer(Genotype) |> strrep(x = 'T', times = _),
         cohort = 'NMOSD Tianjin (n=89)', sex, .keep = 'none') |>
  bind_rows(rs188_1kgp) |>
  mutate(genotype = str_count(genotype, 'T') |>
           case_match(0 ~ 'CC', 1 ~ 'CT', .default = 'TT'))

rs188_2coh |>
  ggplot(aes(y = cohort, fill = genotype)) +
  geom_bar(position = 'fill') +
  theme_pubr() +
  scale_fill_manual(values = c('grey','royalblue','tomato')) +
  labs(x = 'Proportion', title = 'rs188721059 genotype: NMOSD vs HC cohort')

ggsave('mission/NMOSD-SNP/figs/rs188721059_barplot.png', width = 8, height = 3)

rs188_2coh <- rs188_2coh |>
  count(genotype, cohort, sex) |>
  mutate(sex = sex == 'M',
         count = n) 

rs188_or <- rs188_2coh |>
  model_OR_forests(ref_group = '1KGP Beijing (n=103)', alt_allele = 'T')

rs188_sexcor <- rs188_2coh |>
  model_OR_forests(ref_group = '1KGP Beijing (n=103)', alt_allele = 'T',
                   covariate = '+sex')

rs188_2or <- list('None' = rs188_or,
                  'Sex' = rs188_sexcor) |>
  bind_rows(.id = 'covariate')

rs188_2or |>
  write_source_csv('rs188721059_odds_ratio', project = 'mission/NMOSD-SNP/')

rs188_2or |>
  ggplot(aes(or_value, model)) +
  geom_point(size = 6) +
  geom_text(aes(label = str_c('p=',signif(p.value, 3))), nudge_y = .35) +
  geom_linerange(aes(xmin = ci_lower, xmax = ci_upper)) +
  geom_vline(xintercept = 1, linetype = 'dashed') +
  scale_x_log10() +
  theme_pubr() +
  facet_wrap(~covariate, labeller = 'label_both', ncol = 1) +
  labs(x = 'Odds ratio (95% CI)',
       title = 'rs188721059 correlation with NMOSD')

# rs150437294 ------------
snp_id <- 'rs150437294'

rs150 <- read_delim('mission/NMOSD-SNP/rs150437294.tsv') |>
  mutate(Genotype = case_match(Genotype, '0？1？' ~ '1',
                               'Unknown' ~ '2',
                               .default = Genotype))

rs150 |>
  ggplot(aes(`465-510`, `533-580`, color = Genotype)) +
  geom_point()

rs150_1kgp <- query_snp_maf(snp_id, type = 'genotypes') |>
  left_join(pop_1kgp) |>
  filter(subpop == 'CHB') |>
  mutate(sex, genotype, cohort = '1KGP Beijing (n=103)', .keep = 'none')

rs150_2coh <- rs150 |>
  inner_join(tjcvr) |>
  mutate(genotype = as.integer(Genotype) |> strrep(x = 'T', times = _),
         cohort = 'NMOSD Tianjin (n=89)', sex, .keep = 'none') |>
  bind_rows(rs150_1kgp) |>
  mutate(genotype = str_count(genotype, 'T') |>
           case_match(0 ~ 'CC', 1 ~ 'CT', .default = 'TT'))

rs150_2coh |>
  ggplot(aes(y = cohort, fill = genotype)) +
  geom_bar(position = 'fill') +
  theme_pubr() +
  scale_fill_manual(values = c('grey','royalblue','tomato')) +
  labs(x = 'Proportion', title = str_glue('{snp_id} genotype: NMOSD vs HC cohort'))

ggsave(str_glue('mission/NMOSD-SNP/figs/{snp_id}_barplot.png'),
       width = 8, height = 3)

rs150_2coh <- rs150_2coh |>
  count(genotype, cohort, sex) |>
  mutate(sex = sex == 'M',
         count = n) 

rs150_or <- rs150_2coh |>
  model_OR_forests(ref_group = '1KGP Beijing (n=103)', alt_allele = 'T')

rs150_sexcor <- rs150_2coh |>
  model_OR_forests(ref_group = '1KGP Beijing (n=103)', alt_allele = 'T',
                   covariate = '+sex')

rs150_2or <- list('None' = rs150_or,
                  'Sex' = rs150_sexcor) |>
  bind_rows(.id = 'covariate')

rs150_2or |>
  write_source_csv(str_glue('{snp_id}_odds_ratio'))

rs150_2or |>
  ggplot(aes(or_value, model)) +
  geom_point(size = 6) +
  geom_text(aes(label = str_c('p=',signif(p.value, 3))), nudge_y = .35) +
  geom_linerange(aes(xmin = ci_lower, xmax = ci_upper)) +
  geom_vline(xintercept = 1, linetype = 'dashed') +
  scale_x_log10() +
  theme_pubr() +
  facet_wrap(~covariate, labeller = 'label_both', ncol = 1) +
  labs(x = 'Odds ratio (95% CI)',
       title = str_glue('{snp_id} correlation with NMOSD'))

# LD plot ------------
igh_blk <- list(tjcvr, rs146, rs150, rs188) |>
  purrr::reduce(inner_join, by = join_by(ID)) |>
  mutate(rs146563847 = Genotype.x, rs150437294 = Genotype.y,
         rs188721059 = Genotype, `IGHG1-G396R`, ID, .keep = 'none')

igh_blk <- igh_blk |>
  mutate(across(-1, \(x)case_when(x %in% c('RR',2,'2') ~ 'Homo',
                                  x %in% c('GR',1,'1') ~ 'Hetero',
                                   .default = 'WT') |> fct_relevel('WT','Hetero'))) 

g1 <- igh_blk |>
  count(`IGHG1-G396R`, rs146563847, .drop = F) |>
  ggplot(aes(rs146563847, `IGHG1-G396R`, label = n, fill = n)) +
  geom_tile(color = 'black') +
  geom_text() +
  theme_pubr(legend = 'none') +
  scale_fill_gradient(low = 'white', high = 'tomato') +
  labs(title = 'IGHG1-G396R linkage with other IGH block SNPs in NMOSD cohort',
       subtitle = "D' = 0.954")

g2 <- igh_blk |>
  count(`IGHG1-G396R`, rs150437294, .drop = F) |>
  ggplot(aes(rs150437294, `IGHG1-G396R`, label = n, fill = n)) +
  geom_tile(color = 'black') +
  geom_text() +
  theme_pubr(legend = 'none') +
  scale_fill_gradient(low = 'white', high = 'tomato') +
  labs(subtitle = "D' = 0.884")

g3 <- igh_blk |>
  count(`IGHG1-G396R`, rs188721059, .drop = F) |>
  ggplot(aes(rs188721059, `IGHG1-G396R`, label = n, fill = n)) +
  geom_tile(color = 'black') +
  geom_text() +
  theme_pubr(legend = 'none') +
  scale_fill_gradient(low = 'white', high = 'tomato') +
  labs(subtitle = "D' = 0.880")

g1 + g2 + g3

publish_pdf('IGH_LD_block_matrix.pdf', width = 200, height = 100)

## LD value ---------
shelf(genetics)

igh_geno <- igh_blk |>
  mutate(across(-1, \(x)case_when(x == 'Homo' ~ 'T/T',
                                x == 'Hetero' ~ 'T/C',
                                .default = 'C/C')))

LD(genotype(igh_geno$`IGHG1-G396R`), genotype(igh_geno$rs188721059))

## combined OR? ---------
igh_int <- igh_blk |>
  mutate(across(-1, \(x)case_when(x == 'WT' ~ 0,
                                  x == 'Homo' ~ 2,
                                  .default = 1)))

nmosd <- tjcvr |>
  dplyr::select(-`IGHG1-G396R`) |>
  inner_join(igh_int) |>
  write_source_csv('tianjin_89nmosd_5snp')

nmosd <- nmosd |> mutate(IGHG1_G396R = `IGHG1-G396R`,
                         ID = NULL, `FCGR2B-I232T` = NULL, .keep = 'unused')

nmosd

igh_1kgp <- c('rs146563847', 'rs150437294', 'rs188721059') |>
  set_names() |>
  map(query_snp_maf, type = 'genotypes') |>
  list_rbind(names_to = 'id')

igh_1kgp <- igh_1kgp |>
  pivot_wider(names_from = id, values_from = genotype) |>
  inner_join(pop_1kgp) |>
  filter(subpop == 'CHB')

int_1kgp <- igh_1kgp |>
  dplyr::select(sample:IGHG1_G396R, sex) |>
  mutate(rs146563847 = str_count(rs146563847, 'C'),
         across(rs150437294:IGHG1_G396R, \(x)str_count(x, 'T')))

int_1kgp |> write_source_csv('1KGP_CHB_IGH_4SNP')

igh_2coh <-
list('1kgp' = int_1kgp, 'nmosd' = nmosd) |>
  list_rbind(names_to = 'cohort') |>
  mutate(sample = NULL, sex = sex == 'M')

igh_sum <- igh_2coh |>
  rowwise() |>
  mutate(igh_count = sum(rs146563847, rs150437294, rs188721059, IGHG1_G396R),
         all_homo = igh_count == 8)

igh_addit <- igh_sum |>
  ungroup() |>
  mutate(disease = cohort == 'nmosd', igh_count = igh_count/4) |>
  summarise(count = n(), .by = c(sex, igh_count, disease)) |>
  list() |>
  map2(c('','+sex'),
       \(x,y)regress_logistic(x, model_name = 'igh_count', covariate = y))

igh_addit <- igh_addit |>
  set_names(c('None','Sex')) |>
  list_rbind(names_to = 'covariate')

igh_homo4 <- igh_sum |>
  mutate(disease = cohort == 'nmosd') |>
  ungroup() |>
  summarise(count = n(), .by = c(sex, all_homo, disease)) |>
  list() |>
  map2(c('','+sex'),
       \(x,y)regress_logistic(x, model_name = 'all_homo', covariate = y))

igh_homo4

igh_hetero4 <- igh_sum |>
  ungroup() |>
  mutate(disease = cohort == 'nmosd',
         all_hetero = igh_count >= 4) |>
  summarise(count = n(), .by = c(sex, all_hetero, disease)) |>
  list() |>
  map2(c('','+sex'),
       \(x,y)regress_logistic(x, model_name = 'all_hetero', covariate = y)) |>
  set_names(c('None','Sex')) |>
  list_rbind(names_to = 'covariate')

igh_hetero4

igh_homo4 |>
  set_names(c('None','Sex')) |>
  list_rbind(names_to = 'covariate') |>
  bind_rows(igh_addit) |>
  bind_rows(igh_hetero4) |>
  mutate(model = case_match(model, 'igh_count' ~ 'Haplotype_Additive',
                            'all_hetero' ~ 'Haplotype_Dominant',
                            .default = 'Haplotype_Recessive')) |>
  ggplot(aes(or_value, model)) +
  geom_point(size = 6) +
  geom_text(aes(label = str_c('p=',signif(p.value, 3))), nudge_y = .35) +
  geom_linerange(aes(xmin = ci_lower, xmax = ci_upper)) +
  geom_vline(xintercept = 1, linetype = 'dashed') +
  scale_x_log10() +
  theme_pubr() +
  facet_wrap(~covariate, labeller = 'label_both', ncol = 1) +
  labs(x = 'Odds ratio (95% CI)',
       title = str_glue('IGH LD block (4 SNP) correlation with NMOSD'))

publish_source_plot('IGH_LD_block_4SNP_OR_forest', width = 150, height = 100)

g2r_or <- read_tsv('mission/NMOSD-SNP/results/IGHG1_G396R_odds_ratio.tsv')

block_or <- read_csv('mission/NMOSD-SNP/results/IGH_LD_block_4SNP_OR_forest.csv')

g2r_or |>
  mutate(model = str_c('IGHG1-G396R_', model)) |>
  bind_rows(block_or) |>
  mutate(type = ifelse(str_starts(model, 'Haplo'), 'IGH LD block Haplotype', 'IGHG1-G396R')) |>
  ggplot(aes(or_value, model, color = type)) +
  geom_point(size = 6) +
  geom_text(aes(label = str_c('p=',signif(p.value, 3))), nudge_y = .35, color = 'black') +
  geom_linerange(aes(xmin = ci_lower, xmax = ci_upper)) +
  geom_vline(xintercept = 1, linetype = 'dashed') +
  scale_x_log10() +
  theme_pubr() +
  facet_wrap(~covariate, labeller = 'label_both', ncol = 2) +
  labs(x = 'Odds ratio (95% CI)',
       title = str_glue('IGH LD block haplotype (4 SNP) correlation with NMOSD'))

# AUC ROC ---------
shelf(pROC)

roc_result <- roc(igh_sum$cohort, as.numeric(igh_sum$igh_count/4),
                  levels = c("1kgp", "nmosd"), ci = T)

roc_result |> glimpse()

plot(roc_result, main = "ROC Curve for GeneX", col = "blue", lwd = 2)

coords(roc_result) |>
  ggplot(aes(1-specificity, sensitivity)) +
  geom_path() +
  annotate(x = 0, y = 0, xend = 1, yend = 1, linetype = 'dashed', geom = 'segment')

ci.coords(roc_result, x = roc_result$thresholds)

auc(roc_result)
